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NOMENCLATURE A 


Indexing Variables i 


i, 3, 0. 

m 


n 


P 

s 

y 

i 


Integers 

Degree of the denominator polynomial of a 
transfer function. 

Degree of the numerator polynomial of a 

transfer function 

number of input variables 

number of output variables 

An integer varying from 1 to s and identifying a 
particular output variable 

An integer varying from 1 to p and identifying a 
particular input variable. 


Variables ; 

■Ap-q - Data matrix for k th output at q^ instant of time. 

a vi „ j * * •« a i - Coefficients of the numerator polynomial of a 

0 7 KJL 9 HI 

transfer function relating kth output and 
1th input variables* 

1 *’* ,1:, kl n “ Coefficients of the denominator polynomial of 

a transfer function r elating kth output and 
1th input variables. 

c H,i =e J- iT ^ ” Change in the kth output variable due to 1th 

input variable alone at time t = iT. 






- z-transform of change in the kth output 
variable due to 1th input variable alone. 

- Total change in the kth output variable duo 
to all the ’ p’ input variables acting 
together at time t = iT. 

- z- transform of total change in the kth output 
variable due to all the ’p’ input variables 
acting together. 

- A (sxl) vector containing 1 s* total output 
variables in their z- transforms as its 
element s* 

- A (qxl) vector containing kth total output 
variables at qth instant of time as its 
element . 

- 1th input variable at time t = iT . 

- z~ transform, of 1th input variable. 

- A (pxl ) vector containing ’p* input variables 
in their z-transform, as its elements. 

- Sampling interval. 

-A^(m+n +1) vector 'containing 

the coefficients a,, . . . , a. and b,. . , . . 

k1,o> ’ KL,m KL,1 7 

^kl n ^° r ^^put anc3 - fori = 

- A £(m + n + 1 )!P x ij vector containing the 
coefficients a kl Q ,..., m and b kl 

b^ ^ for kth output and for 1 = 
at qth iteration. 



xi 




A 



A {(m + n +1 )P x 1j Sector which best 
approximates the vector X^.* 

A [(m + n +1 (Pxl j vector which best 
approximates t e vector X k at qth instanc 
of time* 
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NOMENCLA TURE B 


•^1 > -^2 



H R 


= cross-sectional area for main steam flow 
and secondary steam flow at the throttle 
valves (in 2 .) 

= constant which relates the flow and pressure 
drop in a turbine with many stages 
= coefficient used in determining steam flow 
through main and secondary throttle valves 
(1b m /sec) 

= heat transfer coefficient of feed water 
heaters (ws/lbl 

= heat transfer coefficient of reheater 
(ws/lb°F) 




h f» b fw ,h fw’ h 0 ~ 


J, k t 
K BEP’ K ap 


enthalpy of main steam at nozzle chest, 
reheater, throttle valve, moisture separator 
and isentropic end points from pressure 
P c and P R , respectively (wsAb) 
enthalpy of saturated water in the reheater 
feed water leaving heater 1 , heater 2 and 
water in the condenser r espectively(wsAb) 
constant to relate mechanical energy and 
vice versa (ws/ftlb,ftlb/ws) 
coefficient which indicates the amount of 
bled steam from the H.P. and L.P. turbines, 
respectively. 
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Q H1* %I2 


9r 

R 

T ni * * 

T R2’ T w 2 ,T w3 


T T 
HP’ LP 

TORQUE 

V , V_ 

C’ R 


v f> v g 


constant in the Callender’ s empirical 

equation relating pressure, density and 

3 

enthalpy of superheated steam (lb^ft / 
in. 2 ws) 5 lb f ft^/in. 2 ) 
pressure in-reactor, nozzle chest and 
reheater (lb/in. ^) 

heat transferred in heater no. 1 and 2, 
respectively (vs) 

heat transfer rate in the reheater (w) 

gas constant (lb- ft'Vcin. 2 lb °F)) 

i m 

time constants associated with heater no. 1, 

heater no. 2, reheaters for flow rate and 

heat transfer rate, and with the H.P. and 

L.P. turbines, respectively (1/sec) 

torque of high pressure and low pressure 

turbines respectively (lb ft) 

f 

net torque at the turbine shafts (lb^ ft) 
total effective volume of the nozzle 
chest and reheater, respectively (ft^) 
specific volume of saturated water and 
steam, respectively, around a pressure 
of 200 p.s.i. (cftAb m ,cftAb m ) 



xiv 


w 1 , w 21 -w^, = steam flow rate at throttle .valve, 

nozzle chest, L.P. turbine, before 
the reheater and after the reheater 
respectively (lb /sec) 

W BHP’ W BLP ~ tied steam frcm high pressure and 

„ low pressure turbines, respectively 
(lb m /sec) 

w FW» w xfP’ W MS = Wa; k er flow rate to reactor, collected 

in heater 2, and separated in the moisture- 
separator (lb /sec) 

w L = a constant wMcb depends on ratio of 

turbine flow required to produce a net 
torque to rated flow of steam 

1 

Wp|^, Wp R = secondary steam flow rate entering and 

leaving the reheater, respectively 

(lb /sec) 
m 

= denotes deviation from the operating point 
= isentrcpic efficiency of the H.P. and 

> r , * 

L.P. turbines, respectively 
~ efficiency correction factors to take 
into account the reduction of rotational 
losses, etc. at reduced loads for H.P. 
and L.P. turbines, respectively 
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Pr, ^2 
Ji 


density of steam in the nozzle chest, 
in the reheater and at the turbine 
outlet (lb/ft^) 

speed of turbine equal to 12o (rad/sec). 



SYNOPSIS 


A method is developed for low order modelling 
of large multivariable systems. The method is applicable 
to both linear as veil as non-linear systems and does 
not require the transfer function or vector differential 
equation of the system to be specified. Onlv measured 
input-output data cf the system at discrete intervals 
of time is required to determine its linear , low order, 
discrete- time model. The method is based cn the use of 
matrix pseudo-inverse to estimate the parameters of the 
model which minimize the sum of the squares of the errors 
between the responses of the original system and the 
low order model at the sampling instants. A recursive 
algorithm is used which reduces the computational and 
data storage requirements and enables successively 
improved model of the syste j to be determined as new 
data from the system is made available. Thus the method 
can be used for on-line applications. 

•An iterative technique is used to convert the 
discrete-time model into a corresponding continuous- 
time model. 

The method is used to obtain low order model of 
a nuclear reactor turbine system. 



CHiPTm I 


1 «1 Intro d uction 

Modern systems are often large, resulting in 
dynamic model which consists of a large number of non- 
linear differential equations* It is common practice to 
linearize the model around an operating point and assume a 
time-invariant behavior* The resulting linear model describ 
the behavior of the system adequately in the neighborhood 
of the operating point and in terms of computational effort 
the system simplifies considerably* In many situations, 
such as a thermal power plant, a nuclear reactor system, 
a process plant or aircraft systems etc., the dynamic 
model is still of such a large order as to be inconvenient 
or impractical for its simulation or control system design. 
Thus, it sometimes becomes important to look for a dynami- 
cally equivalent model which is of reduced order. The 
reduced order model should describe the behavior of the 
system almost as well as more complicated higher order 
model and retain the most important parameters of the 
ori gin al system. 

In the last decade or so considerable attention 
has been devoted to the problem of deriving reduced- order 
models for complex systems. The problem is clearly of 
great interest to those working in the applications area 

£ ^ » \J*~\ 9 
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1 .2 Uses of reduced order models 

deduced- order models are obtained for the 
following reasons : 


(1 ) To simplify understanding of the system, whether the 
problem is analysis, synthesis, or identification and to 
achieve easy simulation of the processes* 


(2) To reduce computational requirements and to make 
better use of the supporting facilities as well as the 
computer itself, in addition to widening the scope of the 
problems which can be handled on a given size of machine. 
Of particular significance to on-line computer control 
applications is the reduction of the computational effort 
of optimal and adaptive controllers by deriving suboptimal 
strategies based on reduced order models^^J, ^5*3* 


(3) To obtain a lower dimensional ’’control law” for 
simplifying the structure of the feedback controller to. 


(4-) To make best use of scanty experimental data, by 
estimating a few parameter with confidence, rather than 
estimating more parameters with less confidence. The 
low order model is then a more reliable predictor tn- 

(5) To generalize results established on a particular 
system to comparable systems, using low order models 
which basically differ in time scale only* This is 
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particularly useful in tolerancing performance for system 

• 

1 »3 Liter at ure s u rvey 

Although most of the work in the area of low 
order modelling of complex systems has been done in the 
last decade or so, the subject has aroused widespread 
interest. In the beginning different investigators 
presented radically new approaches to the same problem. 
Later on several variations and combinations of these 
approaches were developed. These different approaches 
along with their modifications are presented below grouped 
under different categories. 

1 .3.1 Modelling from loop data 

The Initial impetus came from Axelby ^ 9"^ and 
Biernson ^1 0j[ who developed methods for rapidly estimating 
system poles from the open loop Bode plot. This approach 
was formalised by Kan Chen t'h, who suggested that 
dipoles might be estimated from system behavior in the 
high gain region, far-off poles from the low gain region, 
and dominant poles from the cross over region. Kan Chen 
method has been 1 exploited by Can f i el d 2"] for multiloop 
systems, and by Towill U133 in terms of a rational design 
procedure. The basic problem with these methods is that 


check-out £ 



if 

if there are many break points in the high and low gain 
regions, the dominant poles, and hence the transient 
response of the low order model, is considerably in 
error. 

1 «V2 Modelling from the s y stem t ran s fer function 

Since the system transfer function contains all 
the data needed to model the system, there has always been 
effort to use the coefficients directly in a low order 
model. Deriving the system transfer function coefficients 
is time consuming if a computer sub-routine is not available, 
so that for high order systems it is not a preferred 
approach for hand computations* Oldenburger DM neglected 
some of the coefficients on an intuitive partitioning 
basis and applied it to develop a reduced order model of 
a turbire governer system. Siljak^l !?1 , in parameter 
plane design methods, simply neglected the higher order 
coefficients aid varied the lower order coefficients by 
adding a compensator to yield the original system response. 
The method due to Chen and Shieh is based on continued 

fraction expan sicn (CFE) of the transfer function of the 
original sj^sbem. The numerator and denominator polynomials 
in the transfer function are arranged in the ascending 
order, in powers of s, and then expanded into the Caner- 
type continued fraction. From a consideration of the 
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final value theorem, it follows that quotients in the 

expansion are in the order of decreasing significance* 

Thus truncating the CFE at a suitable stage gives the 

reduced transfer function of the desired accuracy. The 

GFE used is equivalent to Taylor's series expansion about 

s = 0 and can be shown to be a special form of Fade' 

approximation for a seal or rational function The 

method ensures that the model gives the correct steady 

i 

state response for a class of inputs, namely^ ^ t , 
i = 0, 1, 2,..*p where 'i' is a positive integer depen- 
ding on the order of the model. However the initial 
transient response may not be good. Shamash^l , 
expanded the numerator and denominator into a Jordan-type 
continued fraction, about s = 0. This method requires 
much less computation than the original method. Ghuang^ $fj-, 
1*0 modified Chen and Shi eh 1 s method. He expanded the 
system transfer function into a Cauer-type continued 
fraction about s = 0 and s = 00 * By alternately expanding 
around s = 0 and s = OO , this method ensures equally 
good steady state and initial transient responses. 

However, biased responses may also be obtained. By 
'fitting 1 more coefficients in the power series expan- 
sion about s = 0 or s = CO , a better stead v state or 
initial transient response respectively, may be obtained. 
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Ghen extended his method to multivariable systems 

by expressing the transfer function matrix as a fraction 
with a matrix polynomial in the numerator and a seal or 
polynomial in the denominator in the ascending order and 
then expanding into matrix continued fraction as in the 
case of scalor systems* Truncating the continued fraction 
at a suitable s tage gives the reduced order multivariable 
model* Chen applied this method to reduce the model of a 
gas-turbine system. 

Vittal Rao and Lamba ^2^3 minimised the integral 
of the weighted error function between a given transfer 
function and the reduced transfer function in the frequency 
domain. This method is of especial use in the non-linear 
systems where the frequency response of the linear part 
of the simplified system should approximate the frequency 
response of the linear part of the original system. 

Shamash l 23 ) extended the continued fraction 
methods for linear continuous- time systems to linear 
discrete- time systems* 


1 .3*3 Eigenvalue approach 

Evans r ^3 is considered the originator of this 
approach as his root locus method brought the determination 
of eigenvalues into the realms of practicability for 
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realistic order systems* The first computerised, state- 
space approach to the problem is credited to Davis cn 
and essentially consists of neglecting eigenvalues of the 
original system farthest from the origin, retaining only 
the dominant time constants. The retention of dominant 
eigenvalues makes the response of the reduced model 
approximate that of bhe original, since the eigenvalues 
neglected make a very insignificant contribution to the 
total response except at the beginning. Relationships 
from the time- solution of the vector differential equation 
of the original model were used to develop a reduced model 
which maintained both the correct proportion of the eigen- 
vectors and the desired eigenvalues. The resulting low 
order model did not preserve steady state values, although 
the initial and middle parts of the responses were good 
representations* Marshall !». by a similar reduction 
procedure and by neglecting the dynamics of fast modes 
of the original high order model, derived a lower order 
model which additionally achieved steady state comparar- 
bility also and offered zero steady state error for cons- 
tant inputs. Following the publication of Davison's work, 
there was much correspondence between Davison and 
Chidambara , but according to Marshall , no new 
method was proposed and the correspondence consisted 
essentially of variations upon the known theme. 
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Chidambara LO accounted for the steady state error in 
Davison’s method by pointing out the terms which cannot 
be neglected as steady state is approached. However, 
Davison 1*283 pointed out that Chid ambar a’ s method while 
giving correct steady state values, did not give good 
transient response agreement. This was because Chidambara 
was essentially finding a reduced forcing function and in 
doing so it did not excite modes of t 1 e approximate model 
in the same proportions as found in the original model. 
Chidambara 1291 •» 1301 countered Davison 1 s arguments by 
giving examples* However, Davison £293 , ^_ 30 ] pointed out 
that the examples given by Chidambara were misleading 
because only one mode was presented In the simplified 
model. He also mentioned that Chidambara’ s method was 
more difficult to use in practice. Following these 
publications, Davison ^3*1 j modified his method by defining 
a diagonal matrix and incorporating it in his original 
method in such a manner that the system gave satisfactory 
dynamic response and correct steady- state response for 
step inputs* Inspite of a number of these articles, it 
was not clear which model should be usd since the examples 
given showed that either method is correct or incorrect 
according to various cases considered. FossardY.323 noted 
this and pointed out that the problem lied in the existence 
of null eigenvalues. He modified Davisons method to take 
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care of null eigenvalues and produced satisfactory steady 
state and dynamic responses for step input. He extended 
his method to general input case also by observing that 
if the slow modes of the system are correctly excited, 
those of the input will be incorrectly transmitted unless 
they are integrated into the system and the reduction 
process made on the augmented system. Newell and Fishei^33j 
applied Davison’s method to reduce the order of an evapora- 
tor model. Seborg and Fisher 5>3 extended the continuous- 
time methods of Davison and Marshall to discrete time 

systems. -An interesting variation to Davison's method is 
that proposed by Kuppurajulu and ELangovan (353 'wherein 
the high order system is replaced by three models, 
successively representing the initial, intermediate and 
final stages of transient response by identifying those 
eigenvalues which are most dominant in the initial, inter- 
mediate and final portions of the system response. 

Mitra’ s ^363 method consisted of orthogonal projection 
onto a subspace of the state space determined by a defined 
matrix. Wilson ^373 determined a low order model by 
minimising a functional of error between the output vector 
of the system and the output vector of the model. The 
algorithm is based on the solution of a Lyapunov 
matrix equation and requires optimization of individual 
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parameters occurring in the particular canonical .form 
chosen* However this method is cumbersome and quite time 
consuming. He extended it to multivariable systems also 
W and used it on a boiler model* 

1 *3*^ Modelling from response data 

This is the general identification problem if we 
probe to find low order models which acceptably describe 
the ‘black box’ under test, through its input-output data. 

In the frequency domain, the historical path commenced 
with weighted least squares method of Levy \39j • He 
fitted the frequency dependent polynomial into the generated 
data by minimising an error function which was a measure 
of error in the fit. The method has been widely used and 
is capable of fitting some quite unusual frequency 
responses ^f0][ , but several difficulities have been noted 
with the method. Sanathana* and Koerner^Vf] found that 
the procedure described does not give a good fit if the 
frequency data spans several decades and have proposed 
instead an iterative procedure. The authors cited an 
example in which the it erativ el y- computed parameters 
differ considerably from those obtained by Levy’s method. 
Sumner has made the same criticism of the Levy method, 

but proposed a technique in which the error between the 
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data and the model is minimised using a non-linear optimi- 
sation method. A further modification made to the method 
by Payne $>3] is to incorporate certain constraints based 
on additional knowledge of the system, such as steady 
state gain. He noted that in the unconstrained system 
there is a tendency for poles to occur on the right half 
of the complex plane and so give unstable responses in 
system known to be stable. He reported, in fact, that in 
fitting several hundred responses about 30 % resulted in 
unstable systems when no constraints were used, but that 
the use of constraints reduced this to only 1 per cent. 
Towill and Payne £,¥+1 studied many systems and came to the 
conclusion that the least squares error criterion did not 
necessarily result in a satisfactory model in many systems. 
However he found this criterion to be most straight- 
forward in terms of application and least computation time 
consuming in most of the cas s* Considerable time may 
therefore be spent in system study and simulation until 
a satisfactory compromise is reached. 

A number of methods minimising some function of 
the error between the system and the low order model at 
sampled response times have been developed. Aider sen 
determined a reduced order model, the response of which 
approaches that of the original in such a manner that the 
mean-square error between the two responses, over a given 
finite internal, is minimised. This is o c<Z5omroli shed hv 
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using the orthogonal projection theorem in the theory of 
linear "vector spaces to minimize the sum of the squares of 
errors between the responses of the actual system and a 
discretized model at the sampling instants* The vector 
dif ferential equation is then obtained from the vector 
difference equation of this discrete model. Sinha 0+6*3 
derived a reduced order model by expressing the model of 
desired order in terms of a pulse transfer function, with 
its coefficients as unin owns* To find out these unknowns, 
he expressed the transfer function in an equivalent 
difference equation form. This equation when written for 
different time instances gave a set of ordinary algebraic 
equations with the coefficients of the reduced transfer 
function as unknowns* The equations were solved for 
these unknowns by calculating Left Minimum Pseudo-Inverse 
which was shown to give the unknown coefficients 
in such a manner that the sum of squares of the error 
between the original and the reduced model responses was 
minimized at the sampling instances* M iterative algo- 
rithm was used and thus the method required much less 
storage of data and much less computation than Andersen’s 
method* Bereznai and Sinha ) +9~\ proposed another method 
based on the minimization of deviation between the 
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responses of the original system and the reduced order 
model in the direction that is everywhere perpendicular 
to the response of the system. He considered various 
error minimization criteria. However, he noted that 
analytical solution to the problem of error minimization 
is available only for the least squares error criterion. 

He suggested the pattern search techniques of Hooke and 
Jeeves ^50]} to solve problems formulated in other error 
minimization criteria. Meier and luenberger l 51 3 
minimised the mean square difference between system output 
and model output when a random process is used to drive 
both systems. Hli £ 523 obtained a low order model from 
impulse response estimates by computing the time moments 
of the impulse response. The moments were then, using 
series expansion theory, equated to the unknown transfer 
function, and the coefficients determined by elimination. 

As the model order increases, considerable difficulty is 
experienced with convergence of the impulse moments for 
realistic test data. Nevertheless, it is a powerful 
technique for use in heavily damped systems. Higgs and 
Edgar (2 53 } developed a general framework for optimum 
model reduction based cn impulse responses of the original 
and reduced single input - single output systems and gave 
example of a heat diffusion system. G-aliana suggested 

a method of approximating linear time invariant 



multi input - multi output systems by a lower order system 
through a weighted least squares minimization of a scalor 
integral functicn of the impulse response matrix over all 
positive time* A set of necessary conditions was derived 
in the form of a number of algebraic matrix equations, the 
solution of which yielded the unknown system parameters. 
Although he gave several examples, he did not apply it 
to any real system* -Applevich' s £553 mei:1:icd on multi- 
variable linear time-invari ant discrete- time system was 
to optimally approximate the given system with respect to 
the sum of the squared output errors of the impulse response. 
He also established existence and uniqueness of the 
solution and in this connection developed certain projec- 
tion theorems. Although he did not apply his method to a 
practical system, he gave an example. 

1 .3.5 Other approaches for low order modelling 

Shi eh and Wei £56 3 P r °P°sed a method which was 
a combination of dominant eigenvalue and continued frac- 
tion approaches. In their method which was applicable to 
multi-input multi- output systems, they expressed the 
transfer function matrix as a product of two polynomials* 

The dominant eigenvalues were used to formulate the common 
denominator polynomial of a reduced order model and the 
numerator dynamics of the reduced order model were 



obtained by matrix- continue! fraction approach* The 
method has the advantage of eigenvalue approach in the 
sense that reduced model is always stable. It also has 
the advantage that the computationally unattractive 
procedure of complicated linear transformation, matrix 
diagcnalizaticn , and steady state value matching, etc. 
involved in eigenvalue approach is replaced by a simple 
technique of continued fraction. Arumug am and Rama Moorty 
'5 7^ used Schwarz canonical form in simplifying dynamics of 
systems described by linear differential equations* The 
method does not require the computation of eigenvalues aid 
eigenvectors* But this method cannot be used for 
multivariable systems. Ravi ehandr an \j)&\ used singular 
perturbation method to reduce the model of a nuclear 
reactor system. 

1 Present work 

In the present work a method has been suggested 
to determine low order model of large multivariable system. 
This method is developed following the guidelines of 
Sinha 1 sf ^63 work. Sinha developed a low order modelling 
technique from input-output data of a linear or non-linear 
system, measured at discrete intervals of time. Low order 
model is determined by using matrix pseudo-inverse 
technique which minimizes the squares of error between 
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the moc 1 el and system responses at the sampling instants* 
However, this method was exclusively for single-input 
single- output systems* Moreover, this method was not 
applied to any real system. Sinha selected a fictitious 
linear system with its parameters chosen in such a 
manner as co simplify the analysis. 

The present method is applicable to multivariable 
systems and can be used for on-line applications* The 
method is applicable to both linear and non-linear systems 
and does not require system transfer function or vector 
differential equation to be specified. Only the input- 
output data of the system measured at discrete intervals 
of time is required to determine the low order discrete 
time model. As in Sinha' s method, matrix pseudo- inverse 
is used to minimize the squares of error between the 
responses of the low order model and the original system 
at sampling instants* 

An iterative algorithm is used which reduces the 
computation and storage of data. The discrete time model 
is converted into corresponding continuous time model by 
application of Smith 1 s L 6 d method • 

The method is used to obtain low order model of 
a nuclear reactor turbine system. 



17 


Chapter U describes the mathematical develop- 
ment of the present method and the algorithm for obtaining 
the lower order model. In Chapter III, the model of the 
nuclear reactor turbine (which has been used here as an 
example) has been discussed and application of the present 
method has been developed. In Chapter IV, the computa- 
tional results considering various orders of reductions 
have been given. Concluding remarks, limitations of the 
application of the method aid scope for further work 
has been discussed in Chapter V. 



CHAPTER II 


P ROBLEM FQRhULATIQK AMD METHO D OF SOL U TIO N 

We oim at the determination of a low order 
model of a multivariable system from its input- output 
data obtained at discrete intervals of time. Such a case 
might arise in modelling of a large system by a low order 
model, which has considerably smaller order dynamics, by 
using directly the measured input -output data. 

Once the algorithm to find lower order model 
from actual input-output data has been incorporated in 
the computer, it will give successively improved lower 
order model as density of the input-output from the plant 
increases. This can be used as an cn-line process. 

2.1 Problem formula tion 

The low order model of a system may be represented 
by a block diagram shown in Fig. 2.1. For the ’p' inputs 
r^ OC) , r 2 C^),*«*r G£), there may be ’S' outputs 

(- 4 ), (*), ... C s (/£) of the model. The relationship 

between the input and output variables is described by 
transfer function matrix H (A), which in effect describes 
the low order model. 

If the system is modelled as low order linear 
model, each input can be treated independently of each 
other. Complete output of the system can then be obtained 





Figure 2.1 Multi-Input Multi-Output System 




by superposition, i.e. , outputs corresponding to' each 
input alone are added together. Thus, the i-th output 
C ± (*) is given by, 

P 

c i (■%) = YL- H i c ? r 3 C4) ; i = 1 (2.1) 

where r j (-%) is the j-th input and H^^is the transfer 
function between i-th output and j-th input with all 
other inputs reduced to zero. For different values of i 
and j, this describes the transfer function matrix with 
each H^(^5)as an element transfer function. Now, the mode- 
lling process reduces to finding coefficients of each 
of the transfer function such that the output variables 
of the low order model approximate those of the original 
system. 

2.1.1 Low order modelling for discrete system data 

Since in practice it is easier to make measure- 
ments at discrete intervals of time than to measure them 
continuously, the system is modelled as low order discrete 
model. The model will be derived so that its output 
approximate those of the original system at the sampling 
instants. Such a model can be described by a digital 
element wherein the output is in the form of pulse 


sequences. 



2*1 


In a digital element, the difference equation 
relating the input and output pulse sequences is, 

e klt 1T l- \l,1 0 KlI (1 - 1)I J- b kl,2 °KLt (:i ' £)T T- •" 

\l,n °KlC (1 - n)I 3 

■ a kl,o r lt li: l + ^1,1 r l\f 1-1 ^ T 1 + * 0 .,% r li ( 1 - £ > T 1 


+• • * + 


where , 


th 


r^ is the 1 input variable* 


th 


C KL "khe ge in the k output variable 


due to 1 input alone • 


'kl,1 »* *•> \l,n &UVA a kl, o’*” * a kl,m 
various coefficients* 


and ai 


etc. are the 


'I is the sampling interval 
i is an integer. 

This linear relation con be interpreted as a 
formula through which the present output number c-^liT) 
can be computed by taking weighted sums of a fixed group 
of input and output numbers. The process of digital 
element is more clearly depicted in the Fig* 2.2* 

The output c kl at any time t = if depends on 
input q. ’applied at t = iT as well as cn the previous outputs 
at times (i~1) T, (i-2)T,... (i-n)T and on the previous 


• * * 
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Inputs applied at times (i~1)T, (i-2)T,. , (i-m)T. The 
relative effects of these variables in the previous time 
values on the present output variable c kl (i T ) is given 
by the weighting coefficients and &kl Q ... 

a KL,m* 

In a linear multi- input multi- output system, 

superposition theorem can be applied; so that the outputs 

corresponding to each input alone are added together to 

give the total output* In other words, 

P 

0* (IT) = Y- “kl an ; * =1,2,... ,8 (2.3) 

1=1 

where , 

(IT) is the k-th total output at t = IT when 
1 P 1 inputs are applied. 

c kl i s k-th output at t = IT correspon- 

ding to 1-th input alone. 

So, in the multivariable, linear discrete system, the 

appropriate difference equation which relates the input 

and output pulse sequences is obtained by taking s um mation 

on both sides of eqn. (2.2) from 1 to p, and is, 

P p 

^LiTD = n 0* CiT] = Y (a ia,o r iL iT l + a KL,i qLu-uil 

1=1 1=1 p 

+ — +a la,m r it Ci - m)T I + I- (b kl,l 

1=1 


+ ** ,+ b kl, n c ia.C (i '* n) ^i ); k = ' 1 >2,...,s 


(2.4) 
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2.1,2 Z-transform analysis of the system 

3h order to describe the dynamics of the discrete 
process represented by the linear difference equations 
given in the previous section, Z- transform analysis is 
used. 


Multiplying both sides of eqn. (2.2) by 6“ iTs , 
where s is the Laplace- transform variable, and taking the 
summation from i = 0 to o® ? we have , 

if s , r , . . v m ~i -iTs 


j? e-“' S - b U|1 o^IO-Dlle 

1=0 


oo 

i=o 


- H b KL,£ OatCi-SJde- 113 
1=0 

t ~ , — i -iTs -tr — 

Y \l,n c klO i “ n) = H 


• • ♦ — 




a, „ . r Xif]e" lTs + 


kl,o r 


i=0 

OO 


^ ®ki,1 r it C 1 - 1 ) T J e + ^T1 ^l^ r 1 tC ; L-2)f}e 


i=0 

'iO 


-ITs 


i=0 


i=0 


DO 

+ - + H a KL,m I l e_iTS < 2 *« 

1=0 

The lower iqd ex in the summation of is zero, 

which indicates the fact that the values of the sequences 
for the negative time are zero. 
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Putting in eqn» (^*5)? Z = e from the defi- 
nition of Z ~ transform and making use of the 'Shifting 
theorem' which states, 


z‘ J r = r \_(i-o)Tl 


and Z~ 3 o^qi] = 


o© 


o© 


£ 'al 11 ! 2 ' 1 - b a,i H z 


-i-1 


i=0 

— • • • — b 


oO 



^> n L- "kl’ 

i=0 

D© 


firlz 


i=0 

~i~n 


OO 

’*,2 

i=0 


-i-2 


o© 


a 


S r 


-f 


•ki.i 21 z 


-i-1 


KUo / 1 

** my <*> 

i=0 

+ a „ r XiT^Z** 1 


-m 


i=0 


a kl,m i — 1' 
i=0 


( 2 . 6 ) 


Factoring out the common summation for each of the terms 
of eqn< (2.6), these result'’, 

00 - -i -1 _2 

"V c kl^^) Z ~ b KL,1 Z , “ b kl , 2 ’ Z ~"" 0 ~ b kl,n Z ^ = 

i=0 


OO 

XL 

r i t ±0 z_i 

(a kl,o 

h *n,i Z_1 + -’- + a ta,m z “ n ’ 

i=0 





Again from the 

definition of Z- transform, 


o (Z) 

KL 

■1 

vUl 2-1 



i=0 




< to 


and 

q < z > 

= XL 




o 

II 



( 2 . 7 ) 
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Using '"hese definitions in eqn« (2.7)? we have. 


c w (Z) (1-b 


kl , 1 


Z -1 - b. 


kl ? 2 


Z 


V,n Z ' n) = 


r i (Z) (a kl, o + ^1,1 z “ + °*‘ +a ki ? ™ z } (2o8) 


m 


or, 


^k! *Hrl 

f^TzT = li (z) = tS* 


o 


+ 


„-1 


*ia,1 2 +, ° ■ +a kl , m 


,-m 


KL,1 l' 1 *..+ a KLi 


7 -n 


(2.9) 


n 


This expression gives the transfer faction between the 
k-th output and the 1-th inputs For different values of 
k and 1 this describes the transfer function matrix H^CZ) 
with each transfer function (Z) as its element* 

To express eqn. (2,9) in the transfer function 
matrix form, wo write it separately for each input variable 
i.e. for each value of index * 1 * , for a given output* The 
total output variable is obtained by using superposition 
theorem . 


P 


P 


"k 


(z) = n °M. (z) = n.\i (z) q (z) 

1=1 1=1 


(2.10) 


Whiting eqn. (2.10) for different values of T k‘ transfer 
function matrix (Z) is obtained. 
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2.1 .3 A particular cas e 

Ai example of a multivariable sysrem having two 
input output variables follows : 


Here 1=1,2 


+ a^ z + * * * + Z m 


So, c k1 (Z) = 


kl, 0 kl ,1 




1 - b Z — b^* p Z — • • • - b. 
kl j 1 ‘ 


IH ( Z) ; k=1 ,2 


,-1 


and c^(Z) 


a K2,o + V,1 Z + -*" +a k2.“ 2 


kl ,n 

-m 


,m 


i-b. p , z' 1 - b z“ 2 -. ..-b„ 

k2,1 


-n 


r 2 (Z) ;k=1 ,2 


”k2,2 ~ k2,n 

These are the output variables due to inputs r^ (Z) and 
r 2 (Z) respectively. The total output variable is, 

C k (Z) = G k1 (Z > + C !2 (Z) 5 k = 1,2 

Thus the system is represented in the following transfer 
fun c tic., matrix form - 


c 11 (Z)+c 12 (Z) 


c 21 (Z)+c 22 (Z) 


C,(Z) 


C 0 (Z) 


1 


1 , 0 ^ 1 1 , 1 Z ~ + *** +a ' 


7-m 


-1 


1 -b. 


H- 


11_ i.2 J o +a 12,1 2 + * ,,+a 12, m 


-n 


^“"^1 2 "i Z" 1 - • • * -b-, 


x.. 

Z-“ n 


1 1» 1Z " — **- b 11 ,n^ 12,1 ^ “••*" u 12,n 

a 2l , o +a 21 , Z" -1 **- * * + b 21 , n Z ®22 , o +a 2 2, 1 Z 1 + * * * ^22 ,m Z 

1 —bo-i ^ Z ^ -<* « » -b 0 h ^ Z 


-m 


or 


21 ,1 * * 21 ,n 

C (Z) = H( Z) R (Z) 


1 _b 22 , 1 Z " 1 ' • * # “ b 22 ,n Z 


-n 


1 

L: 

** ^ 
:.ri(2) 


r 2 (Z) 
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ZsXjb- -_G.ener.al transfer function m atrix 

As in the section (2.1 o 3), the eqn. (2-9) 
may bo used to write the transfer function matrix cf the 

•tM" 

model having ? p* inputs and 's’ outputs. So, the reduced 
order, linear discrete system may be expressed in the 
pulse transfer function matrix form i 

C (Z) = H(Z) R(Z), (2.11 ) 

where, 

C T (Z) =[C 1 C 2 ... (2.12) 

R T (Z) r 2 ... r^, (2.13) 


H(Z) = 


U off. . .. W /-' 1 i w 

Q-a.ip* ‘ ”>irCL'i.x,rvrfc' "* tX-u'P.a* • '*■ ~V- a -2.^Yyi t: £ V " t ' 

>-k 2 i ar^i -4 _ v-h i -l •*-* l ^ 

' 

l * 1 r 


r** 


***,*** 

(2o1^ 


of from eqn. (2.4-), in the form of the equivalent 
difference equation ; 
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°k,i 


Pm p n 

^ b kl,3 C kl,i-o 


k = 1 ,2, 0 o e S (2.15) 


n p 


0=0 1=1 


\ 

a n n . r + / / b, . c, . 

^->3 1,1-3 £lj KL,J KL,i-3 


k = 1 ,2, . • • , s (2.1 6) 


Letting i range from zero to some integer ' q f , 
eqn. (2.16) may be represented by the matrix equation : 


V ** = %, q 


( 2 . 17 ) 


where : 


V _ C (a k1,0 \2,0 •" ( Vl,1 a k2,1-” a kp,1 ) -” Ca i 


kl ,m 


i, _ ... a, ) (b. j. r«b. „ ^ ) (b „ b _ *. *b, 0 ) ... 

k2,m kp,m' kl ,1 kP,1 ' kl ,2 k2 2 kp,2' 


(b . b. _ .... b, )1 

kl ,n k2,n’ kp,n -* 


(2.18) 


% “ ft ,* %£'***& 


(2.19) 
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( 2 . 20 ) 
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Vector is termed ’Parameter* vector and 
matrix A^ is termed ’Data’ matrix. 

The problem, then, is the determination of 
(m +n + 1 )p components of the Parameter vector X^. from 
the matrix eqn. ( 2 . 17 ) the left hand side of which contains 
the Data matrix, obtained from the input- output data at 
sampling instants. Its right hand side is the vector of 
total output variables, all of which are known. Therefore, 
the problem simply appears to be that of the solution of 
a set of linear algebraic equations. However, in practice, 
the number of rows in the Data matrix is much larger than 
the number of columns. This is particularly so when the 
data is taken over a large interval and the sampling rate 
is high. Moreover since we want the order of the reduced 
model to be as small as possible, the number of columns 
in the Data matrix will naturally be small. Because of 
this reason we get mere number of equations than are the 
number of unknowns, and the set of linear equations repre- 
sented by ( 2 . 17 ) does not possess a unique solution, unless 
m and n are the actual orders of the polynomials of the 


pulse transfer functions. As this is not the case, by 

hypothesis, we can determine that parameter vector, 

■/% 

X^ , which gives the ’best' fit. One suitable criterion 
for an optimum fit is the minimization of the sum of the 
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squares of the errors between the responses of the actual 
system and the reduced model at the sampling instants* 

This is the criterion considered in the present work. 

2.2 Method of solu t ion 

As stated earlier, the problem does not have a 
unique solution. In other words, since the Data matrix 
is rectangular, it does not have an inverse in the ordinary 
sense. However, for 3 > (m+n +1) p, the least-squares 
solution to eqn. (2.17) .is given by, 

A LM 

X, = A C (2.21) 

“k kq ~q 

LM 

where is called 'Left Minimum Pseudo Inverse* of the 

matrix and is defined by, 

LM r LM T -1 

v = [ ( v } v] • ^ (2 - 22) 

The properties of the matrix pseudo-inverse have 
been discussed thoroughly by Penrose L^7"^ and GrevilleQ+81 . 

A description of matrix pseudo- inverse is also given in 

<\ 

the appendix. There it has been shown that X^, from eqn. 
(2.21) is determined in such a manner that the sum of the 
squares of error between the right hand side of eqn .(2.1 7 -) 
(which is the measured response of the actual system) and 
its left hand side (which is response of the low order 
model) is minimized. In other words, pulse transfer 
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functions in the transfer function matrix (2.14), with 
the co"fficients determine^ from (2.13), y”* eld respon- 
ses of the low order model which approximate those of the 
original system. 

2.3 Algorithms based on the matrix pseudo- in verse 

In practice, the Data matrix is quite large and 

LM 

it is not feasible to calculate directly. This is 
particularly the case when a large amount of data must be 
stored in order to minimize the mean-square error over 
a sufficiently large interval, especially when the samp- 
ling rate is high. 'This difficulty can be overcome by 
using the following recursive algorithm, originally 
proposed by Albert and Sittler [.59*1 and modified by 
Wells £ 60 ] f 0 p the case where a row is added to A 
every time, as happens with each pair of additional input 
output data: 




3 ^ 


and let 


'k,q+1 


(2.2?) 


“k,q.+1 


The recursive algorithm that follows is divided in two 
phases* The first phase is the starting phase* !a£ 
start from q = 0 and go upto q < (m + n +1) p and 
use the following recursive relations ; 


-K 1 +1 ^k£L 


V \ A q+i (Ck ’ cl+1 " ^,<1+1 “k,q } 


a Q, a 

"It ,q+1 k,q “k,q+1 


(2.26) 


where 


:,q %,q+lC Q k,q -Sfc.q+l] 


(2.27) 


■%t,q+1 %q %t 5 q +1 


and P k,q+1 = V 


T 

Pk,q ^c,q+l][ P k,q % r q+l} + [\,q%,q+l] 


-k,q+l Q k ? q ■ § k,q+l r 

L k ? q* 


%t,q+1*] 


k, q — k q +• 


-l1l Q k,q ~k,q+ll t 1Ha k,q-fi P k,q %,q+ l] 


r I 

C \,q+1 ^k,q \o.+1 


(2.28) 


A 

\fith the initial conditions X = 0, (1 = I and P = 0. 

— o “ o ° 

The second phase follows the first and starts for q>(m-fo+1)p. 
The recursive relations in this phase are ; A 

m % 


Y _ y _ k.q-^c„q+1 ^k .q^1~ -g 

k q -H ^“k ,q i+a d P a 

1 + -%c v q+1 F k,q %t,q+1 


Xi .SL„t! *1 ~‘kq) 


(2.29) 
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where 

P = P 

k,q+l k,q 


[ P k,q -^jq+^kjq %,q+l] 
+ %c 5 q+1 P k,q %,q+1 


(2.30) 


2 «U- Do to rmin ati. on of . the equivalent cent In nous- time model. 

The algorithm given previously determine the 
low-order pulse transfer function matrix H(z) for the 
system. The problem now is to determine the corresponding 
continuous-time system transfer function matrix H(^). It 
should be noted that H(s) is not the Laplace transform 
of the impulse response which has the z-transform H(z). 

The relationship between H(s) and H(z) depends upen 
the farm of the input between the sampling instants* 
Andersai(V?j| has assumed that the input is held constant 
between the sampling instants. This is equivalent to 
the introduction of a zero- order hold, so that H(z) may 
be regarded as the z-transfer function of H(s) preceded 
by a sampler and a zero-order hold. If the input to 
the system is a unit step this does not introduce any 
errors at all. 

Smith £ 61J has given a method for estimating 
the Laplace- transform coefficients of a dynamic system 
from sampled-.! at a values of its input and output. 

Advantage of this method is that it is quite straight 
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forward and can be easily programmed in a computer as 
the algorithm is iterative in nature <> The method is 
applicable to the systems with noise. Its application, 
on nod. so- free systems is describ.d in section 2.4-.1 . 


2.4-.1 . ... Description of Smith 1 s m ethod 

In this method the dynamic system is assumed 
to he linear with a Laplace transform of the form, 


H(s) = 


04 


+ 


«<-] s + 


+ cC p s 


P 


£ 0 + S + . .0 + (3 s a - 


(2 .1+7) 


The coefficients 06 Q5 CX p and p 0 ,...,$ q 

in this transfer function are not in own. The pulse 
transfer function H(z) which exactly relates the input 
and output sampled data sequences is 


A 

H(z) 


A 


a n + a- z' 


+ au z 


-m 


A A A _n 

b + b, z + ... + b z 
o1 a 


(2.4-8) 


In general each of the coefficient in this 


pulse transfer function is a function of coefficients 
cf the transfer function (3.4-7) • In other words 3 
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/s A 

a o = a o ^ o’ 1 »• ** » <X p - P 0 > (3 q ) 


^ A 


m m 

A A 


a m ( ^ o’ *** 1 ’ * ‘ * ** p » P o s Pi’ 0 * 0 ’ P q> 


b o b o ^ ^o’^l »***»^>< p» Po’ Pi’***’ P cP 


( 2 .^ 9 ) 


A A 

b n = b n 0 » ^*<1 J *•*» p » Po 5 Pi’***’ Pq) 

Now, tho pulse transfer function (2.4-g) must be 
of tho same form as the pulse transfer function H(z), 
estimated from the data. If it is not, then the form 
of H(s) is changed in such a manner that the forms of 

/V 

the nulse transfer functions H(z) and Hf z) are exactly 
similar. 


The pulse transfer function H(z) ustimated 


from data is given by, 


H(z ) 


3-q ~ >r Q-i z 1 . . + a^ z 

b Q + z~^ + ...+b n z“ n 


(2.50) 


The coefficients, a Q> , . . » , a and b Q , b_j,..., 
b n of this transfer function are known® 

Now equating the coefficients of H(z) and H(z) 
for corresponding powers of Z, we get a system of 
equations, withc * 0>c *1 »* •• an ^ Pq’ Pi’***’ Pq 35 
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unkn own s , 

a o |3 q ) = a o 


r 

0 ,0<: 1 ’* ** »^pj P 

o’ Pi »••• 5 P 

•Q^ 

11 

A 

Q J oC 1 , C< p , 



b 

0 

• 

P o ’ P 1 » • • • » 

= b 0 

• 

• 

A 




b n 

(o< 0 , o< -]>•••> oC. p ? 

p o’ — > 

<V = b „ 


Solving eqns. (2.51) fcr the unknowns, o(. , Q<. , 

o 1 

oCp and. |3 o , (3^,..., {3 , the coefficients of the 
transfer function H(s) are obtained. An example of 
this method is given in section 2. 4.2. 

2 .4- . 2 n example of the am Ti cation of Smi^h* s metho d 

Suppose that the pulse transfer function 
obtained from sampled data is, 
b 0 +b 1 z" 1 +b 2 z“ 2 

H(z) = pj P (2.52) 

1 + a^ z + ap z~ 

Here, the values of coefficients b Q? b^ , b^, a^ and 
ag are known* The form of this transfer function 
suggests that H(s) to be estimated from it should have 
the form m , 
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oLo + C?C 1 S 



( 2 . 53 ) 


< 2 . 5*0 



( 2 . 55 ) 

( 2 . 56 ) 

( 2 . 57 ) 



As can be seen, there are four unto own s , ^ , 

c/.p and and five equations relating these 

variables are available. Thus exact values of all of 
these coefficients cannot be determined. To get approximate 



b-0 


values, 
Step I 
Step II 

Step III 

Step IV 

Step V . 


of following iterative scheme may be used, 

- Find^ from eqn. ( 2 . 58 ) 

- Assume a value of p 2 and f ind^4 Q from 
eqn. (2.^0, 

- Combine eqns. (2.56) and 2.57) to eliminate 
trign ometrical expressions and by substituting 
the values of , {5 2 and ^ 0 , findc<^ , 

- Substitute values of < 3 ^, 0 p 2 in eqn. (2.55) 

and get a value ofoc • 

o 

- Compare the values of^ Q obtained in Step II 
and step IV. If they are close enough, stop, 
otherwise, assume a new value of (3 and 
repeat the above steps. 

A computer program can be made to implement 


the above scheme 



CHAPTER. Ill 


APPLICATION TO A NUCLEAR REACT OH TURBINE MODEL 

3jl 1 Br;Lof_A e. script ion of th.e_nncl.par re act or turbine model 

Tlie development of the model of the nuclear 
reactor turbine given here is based on IBM report £62“\ and 
reference 

The block diagram representation of the turbine 

* 

coupled to the nuclear reactor with its various elements 
is shown in Tig. 3.1 • The steam turbine assembly can be 
divided into various subsystems like high pressure, 
intermediate pressure, low pressure units, feed water 
heaters, etc. The well known process representation for 
this part of the modified Rankine cycle is represented in 
Fig. 3*2. The steam flowing from the reactor through 
the main steam valve is just dry, hence the steam exiting 
from the high pressure turbine will be wet, as indicated 
by the condition of steam at point 2 of Fig. 3*2. The 
wafer part of the wet steam is removed in the moisture 
separator and the dry steam at pressure p^ is superheated 
to temperature % in the reheater by using secondary steam 
drawn from the main supply line through a Secondary Steam 
Valve. Tliis superheated steam expands through interme- 
diate pressure and low pressure turbines and then enter*-* 
the condenser. The condensed steam, i«e. feed water, 



flows through the feed water heaters before it enters the 
steam C rum* The feed pump 4 n between the two heaters 
raises the pressure from the condenser pressure to the 
reactor drum pressure. For the above system the model 
and plant dynamics are described below. 

3*2 Model an d tho Plant dynamics 

The model was developed by the application of 
well known physical laws and a few pertinent empirical 
relati on ships* Various assumptions made in the develop- 
ment of the model are presented first, followed by the 
model itself. 

(1 ) The relationship between enthalpy, specific 
volume and pressure (i.e. the equation of state) for the 
superheated steam follows Callender's empirical relation- 
ship for all non-constant pressure processes. 

(2) The I.P. and L.P. turbine dynamics are 
lumped together, as both have a correspondingly slower 
steam flow dynamics as compared to the H.P~ turbine. 

(3) The flow through the throttle valve is 
assumed to depend cn the square root of the upstream 
pressure and density of the steam and linearly on the 
area of flow. 

(*+) The regenerative bleed flow depends linearly 
on flow rate. The flow-time relationship for the turbines 
is represented by a first order lag. 
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(5) In order to take into account the reduction 

of rotational losses? root and tip clearance losses for 

the TI.P. and L.P. turbines at reduced flow rates, 

Vq and H , are multiplied by correction factors 
‘HP Lp 

ft and ^ , respectively, which are assumed to 

t CPH ? 1 CFlP 

depend linearly on the flow rate. 

(6) The bleed flow is tapped right at the end 
of the high pressure turbine, and thus the entire flow 
through the turbjne participates in producing torque 
while only one-half of the bleed flow is assumed to 
produce torque in the L.P. turbine. The second assumption 
is due to the fact that bleed flow is removed at various 
points in the I.P. and L.P. turbines. 

(7) Heat exchange in the reheater is assumed to 
bo perfect and the dynamics of mass balance and energy 
balance are lumped at a single point. Pressure drop in 
the rchcater is assumed to be zero. 

Hith these assumptions the system equations 
representing the turbine dynamics are given below. The 
variable names are explained in Nomenclature B. 

(i) Nozzle Chest 


w 1 h s ~ w 2 h c ^ J p c d .fc, 1 


? V 
c c 


772 dt~ 


1- k J 
i 


( 3 . 1 ) 



d fc = W 1 - W 2 

dt v 

c 


w 


1 


c* A p 
1 1 


w,. 


A* <?„ P„ - P H f 2 ) 


% = (h 2 - bjO/h^ 

p = 1 

^2 x" + (1-x) V|« 

v f = 0 • 01 3*+ 

v g = 2.288 - 0.0166 (p R - 200) 

h = 12.57 X 10^ - 27 (p - 1000) 

S 

p o =/ C (k 1 h c - V 


(ii) High pressure turbine 

dw p 

T w 2 dF = (V 2 ' W BH > ' W 2 


W„ 


BHP 


BHP W 2 


(iii) Moisture separator 


W MS = 


! 



<W 2 ' K BHP V ' W 2 
h 2 - h f ii 

— ‘wl 

hfg 2 

355 + 0.M+ (P R - 200) 


x i .055 x i o 


5 
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(3.2) 

(3.3) 
(3 At-) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 

( 3 . 10 ) 

(3.11) 


(3-1 2) 
( 3 .1 3 ) 


(3-14) 
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h fcr = W - 0.4- (^ - 200) x i .055 x io 


5 


(3.15) 


(iv) Re heater 

(a) Main steam 


R 




dt 

d hR 

dt 


Tr 


4 * 


+ w 2 h g - v 3 


P r 

J R R 

P R = ^ R ^ k 1 h R ~ k 2^ 
1 /2 


w 0 = k„ 


R^ 


Jp R (W 2 “ 




R R 


(3.16) 


^ 1-KjJ ^ 3#1 


( 3 . 18 ) 

(3.19) 


h = [1198.4 + 0.c4 (p^ - 200)“} X 1 .055 X 10 3 ( 3 . 20 ) 


(b) Reheater steam 


_ WpR - W 1 R 

^ 


(3.21 ) 

d Or _ »p R + Vi 

dt 2 1 if (T s - %) 

d X R2 Tl 

%_ 

T R2 

(3.22) 

W PR = C 2 ^2 p 


(3.23) 

T g = 544.6 + 0.125 (p - 1000) 

+ 460 

(3.24) 

T R = Pr^Pr R 


(3 .25) 

R = 201,7/0.365 ( 520 + 460) 


(3.26) 
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(V) I 

*P . and L. P 

. turbine 


dt 

: ^ ” k blp^ 

t 

w 3 

3 r w3 

(3-27) 

W BLP 

= K BLP w 3 


(3.28) 

(vi) 

Heater 1 



dt 

_ 9si 

t hi t ^ 

h - hi 
+ 0 IW 

T 

HI 

(3.29) 

Q H i * 

He (w* + 

Fw v hp 

\p> 

(3. 30) 

(vii) 

Heater 2 



dh fw . 

_ %.z 

f ^fw " h fw 

(3.31) 

dt 

%2 w Fw 

T H2 

cO 

II 

H PV ^ W ELP 

fw ie + w Jr ) 

(3.32) 

dw HP 
ST : 

t 

w 2 **■ *w 2 

4 . W PR ' W HP 

(3.33) 

%2p 

T H2p 


(Tdrii) kvork output equations 


'HP = W 2 (h = - h 2 5 

K i 

(3.34) 

LP ~ Hlp W 3 hR ” ^ 

(3 .35) 


h' 2 =£l067 + 0.37 (P^ - 200) - O.OOH (p^- 200) 2 - . 

O 

o«l (p - iooo)Jx 1 .055 x 10 (3.36) 

bj«. = 898 X 1 .055 x 10 3 
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* 


^ HP ~^GPHP ^ HP 

(3*37) 

LP CFLP *4 LP 

(3*38) 

. * . *■ 
n*=n = 0 .® 

LP HP 

_ w£/w* - w L _ 

CFHP (v'/ w *) ) 

w'/w* - W L 

l cflp = ; . , * , 

(w^/w^) (1- w ) 

(3.39) 

(3.4-0) 

K t = 778/(1 ,075 X 1 O 3 ) 

(3.4-1 ) 

TORQUE = T + T 

HP LP 

(3.4-2) 

TORQUE X fi 

POWER = 

500 

(3*4-3) 


A closer look at the equations reveal that 
equations (3.1), (3.2), (3.10), (3.16), (3.17), (3-21), 
(3*22), (3.27), (3.29), (3«°1) and (3*33) are differential 
equations with the rest being algebraic relationships. 
Thus, in the turbine model there are 11 nonlinear diffe- 
rential equations and 32 algebraic relationships. 

They can be represented in the state variable 

form as, 

i - f(X, U) (3.1*) 

where, 



T 

statevector , X = (X^ , Xg ... 

= (h c M« f Q i> H Ifc h R % l4 w h fw W RP 

( 3 .^ 5 ) 

input vector, U T = (1^ U 2 ) 

= (A, (3A6) 

3.3 Low order model of the plant system 

The method developed in the present work w as 
applied to obtain lav/ order model for the nuclear reactor 
turbine system described above* Although,,, for this system 
the model equations were available, in actual practice, 
there is no necessity for these equations to be present 
in order to obtain low order model for the system. Measure- 
ment of system responses at discrete intervals of time due 
to known excitations to the inputs is sufficient to deter- 
mine low order model for th* system. In the present work, 
this data was generated from the model equations. 

As can be seen from the equations describing the 
plant dynamics of the system, it is not easy to find out 
the exact order of a transfer function relating an output 
variable to a particular input variable. However, looking 
at the number of differential equations in the svstem, 
which is eleven, it can be assumed that no element in the 
transfer function matrix can be of order more than eleven. 
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It was decided to model each element in the transfer 
function matrix as of second order and go for higher 
orders if the responses of the model did not match with 
the responses of the original system. However, strategy 
was to go for even lower order model and find out the 
lowest order model possible. 

In practical situations if plant system equa- 
tions are not available, one can decide on the order of 
the transfer functions and adopt the above procedure to 
obtain low order model. 

Having decided to model each element in the 


transfer function matrix (2.14) as of second order, the 
most general second order transfer function with m = 2, 
h = 2 in eqn. (2*9) was selected. Thus a typical 


second order transfer function will look like, 

a ( , _ a kl,o + a kl,1 2-1 +a kl,2 2 2 

H kl (z) - r—T —£\ ^2 

1 “ b kl,« 2 " \l,2 Z 


A 3*3 percent step change in the main throttle 
valve area (Aj) and a 10 per cent step change in by-pass 
valve area (J^) were made for response studies. Changes 
over the steady s tat e valves in the output variables due 
to both of these inputs acting together, as well as 
individually, were taken with the sampling interval 
T = 0*1 sec. and using 200 samples* 
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These output values were put in the Data matrix 
(2.20) and equation (2.21) was solved for unto own 
coefficients by using the algorithms given in section 
2.3* The computations were carried out on DEC 1090 
Computer System. 

The coefficients of the pulse transfer functions 
were converted into corresponding continuous transfer 
function coefficients by the method described in section 

2 «U- . 1 » 

Other low order models for the system with 
transfer function matrix having all first order transfer 
functions as its elements, and with transfer function 
matrix carrying a combination of first and second order 
transfer functions were also derived. Discussion of 
responses of these various models is given in Chapter IV. 







CHiPTm IV 


DISCUSSION ON THE RESULTS JflD SELECTION OF 
LOW ORDER MODEL 


Intro d uctio n 

The present method was applied to obtain low 
order models of the nuclear reactor turbine system 
described in Chapter m. In t^e present chapter 
different low order models for the system are presented. 
Responses have been obtained for - 

(1 ) Reheater temperature due to 3*3$ step change in 
the main throttle valve area and 10$ step change in by- 
pass valve area, 

(2) Low pressure turbine torque due to 3*3 % step change 
in the main throttle valve area and 10 % step change in 
by-pass valve area. 

(3) High pressure turbine torque due to 3*3% step change 
in the main throttle valve area and 10 % step change in 
by-pass valve area. 

0+) Steam flow rate due to 3*3$ step change in the main 
throttle valve area. 

3 #3$ and 10$ step inputs are chosen only for the 
the sake of comparison of responses cf the low order model 
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and the actual system, the responses of which are known 
for these inputs 1623 - The method at such does not 
put any restriction cn the input. 

4.1.1 Reheater temperatur e 

The four models obtained for this response 

are - 

Model I - Both main and by-pass throttle valveshaving 
second order transfer functions. 

Model II - Main throttle and by-pass valveshaving II and 

I order transfer functions respectively. 

Model III- Main throttle and by-pass valveshaving I and 

II order transfer functions respectively. 
Model IV - Both main throttle by-pass valveshaving I 

order transfer functions. 

These models are presented in table 4.1 (a) 
and their responses are shown in Fig. 4.1 (a). In this 
figure, it can be seen that the best overall response is 
obtained by model I. However, the best peak response 
is obtained by model II. This is because of the ncn- 
oscillatcry nature of the first order transfer function 
(used to model the behaviour of by-pass valve) super- 
imposition of which tends to reduce the overall peak 
overshoot. The transient responses of Models III and. IV 
are considerably in error. Magnitude of steady state 
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error (.07%) is same in all the models. 

Tabic 4*3 lives time constants (I order irons, 
func.) '"nd settling time (II order trais* func.). Table 
4.4- 1 ‘ivjs the GPU time for model generation. 

Computer Processing Unit (CPU) time for the 
various models are shown in Table 4.4. 

4«1 .2 Low pressure turbine torque 

The models obtained for this response are : 

Model I - Both main throttle and by-pass valves having 
second order transfer functions. 

Model II - Transfer functions of main throttle and 
by-pass valves of second and first order 
respectively. 

Model III- Transfer functions of main throttle and by-pass 
valves of first and second order respectively. 
Model IV - Both main throttle and by-pass valveshaving I 
order transfer functions. 

These models are presented in Table 4.1(b) 
and their responses are shown in Fig. 4.1(b). As can be 
seen from this figure* the best transient and steady- 
state (.07% error) responses are obtained by Model I. 
Response of Model II gives much higher steady-state 
error (.21%). The transient response is also not as 
good# Much poorer transient responses are obtained by 
Models III and IV. However, the latter gives much 
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smaller steady-state error (.2$ in model IV compared 

to .42^ in Model III), lime constants, settling time and GPU time 
are v.ven in Table 4.3 and 4*4 respectively. 

4. 1_* JL High pressure turbine torque 

Four models are obtained for this response and 
they are ; 

Model I - Both main throttle and by-pass valveshaving 
second order transfer functions. 

Model II - Transfer functions of main throttle and 
by-pass valves of second and first order 
respectively. 

Model III- Transfer functions of main throttle and 
by-pass valves of first and second order 
respectively. 

Model IV - Both main throttle and by-pass valveshaving I 
order transfer functions. 

These models are given in Table 4.1 (c) and 
their responses are shown in Fig. 4.1(c). The figure 
shows that the best steady-state as well as transient 
responses are obtained by Model I. But the best peak 
response is given by Model II. This is because the 
peak overshoot is reduced by the ncn-oscillatcry behaviour 
of the first order transfer function which is used to 
model by-pass valve characteristics. Model III gives 
higher transient and steady-state (.6 <$) errors, than 
the other two models. However, the highest steady-state 
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and transient response errors are given by Model IV. 

Time constants, settling time and CPU time are given 
in Table 4*3 and 4 .4- respectively* 

4.1 , b- Steam flow rate through nozzle chest 

Two models are considered for tbis response ; 
Model I - Second order transfer function 
Model II - First order transfer function. 

These models are given in Table 4.2 and their 
responses are shown in Fig. 4.2. As can be seen from 
tho figiire, Model I gives the best transient and steady- 
state responses. Model II gives much more transient error, 
however, steady-state responses of these two models 
do not differ much (.10% error in Model I and .15% 
error in Model II). 

Table 4.3 gives time constants' (I oider -trans* 
fume.) and settling tirjg time (II order trars* func.). T 
4J+ gives the CPU tine _ for model,, generation. 

4 a ,2 Selection of a 1 ow or d er model 

In the previous section several low order 
models for the nuclear reactor turbine system were 
obtained. In t 1 is section suggestions have been made 
regarding the selection of a particular model. 

The low order model is selected on the basis 
of application which it is expected to he put to. In 
the on-line applications like adaptive control of a 
system, the need is to generate in smallest possible 



time a 1 ow order mcdeH, which will describe the responses 
of the original system with sufficient accuracy* Here, 
the tjjue required to generate the low- order model is 
the most importaii; factor. Thus a model with smaller 
number of parameters will be given preference over a 
model with more parameters even .though its response may 
not be as accurate as that of the more complex; model. 
However, in the situations where the need is just to 
obtain a low order model on the system, and to continuously 
refine it as more information on the system becomes 
available in the form of input-output data, the time 
required to generate the model is not a very important 
factor. In these cases, the low order model which 
gives the most accurate responses will be selected. 

Based on these considerations we will now select 
low ordji- models of the nuclear reactor turbine system. 

4 .2 .1 Most appropriate mode l for on-line ,, an.pj.ic atj_q?ja. 

As can be seen from Table 4*3, generation time 
of Model II for all the four responses - reheater tempe- 
rature, low pressure turbine torque, high pressure 
turbine torque and steam flow rate through nozzle chest, 
is around 6C$ of that required to generate Model I. 

Thus for all the responses, Model II is selected and 
the low order model of the system for these responses 
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is, 


1 ? 

I 


"3*9896 (1 + 0.9681 s) 

80 21 12 n 


s 2 JO. 6624 s + 0.230’'- 

1506.I1O: (1 + 0.6605 

s + c ,0834 

s) 99-3512 


s 2 + 0.6624 s JO. 23 ( 4 - 

s + O.ORi 8 

1 

■ c 1 


2511 . 87 $- (1 + 0.6831 s) 

1 0.03+ 

3 i 

r . 


s 2 +0.6562 s + 1 .2061 

5*3287 

s+O.OB+O 

s +0.0830 ' 

O 


"2 


-4 


where, 

°i 

% 




Reheater temperature over its steady 
state value 

Low pressure turbine torque over its 
steady state value 

Ugh pressure turbine tcrqu'"' over its 
steady state value 

Steam flow rate through the nozzle chest 
over its steady state value 
Main throttle valve area 
By-pass valve area. 
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4.2.2 Host appropriate model far off-lin e application' s. 

Far off-line appli aticns, time r^iuired to 
generate a model is not the most important factor. The 
model of reasonable order which best describes the 
system responses are selected. So, for the nuclear 
reactor turbine system, Model I is selected for all 
the responses. Hence the most suitable model for 
this system is, 


1 

C 1 


~9*02250 + 0. 8614s) 

8.7032(1 40.1297s) 



s 2 + 0.6952s +0.2959 

1F+ 0.6830s-4)o2314 


C 2 


1 521 . 738 I (1 +0.44 32s) 
s 2 40. 6222 s + 0.26 01 

1 00.6261 (1 40.41 18s) ; 

s 2 -6.6248s -KJ.1 936 ■ 

>r 



2499.0842(1 +0. 5369s) 

11 .5254(1 +0.671 3s) 


°3 


s 2 -O. 672 I s + 1 .2544 

s 2 -io.63*-s+0.8ft-9 

> k 2 

L 

C 4 
. - 


5*41 *4(1 + 0 . 2967 s) 
s 2 -f0. 6971s ^0.7225 

u 

O j 

i 




IV 


8. 979*+ 
s -+0*081 9 


8.1 312 
s 40.0806 
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Table 4.1 (b) 

Low Pr e ssure Turbi n e Tor a ue 


Input i 

able j 

Model \s l 

Change in the ma'.n i Change in the by-pass 

throttle valve area i valve area 

1 

■) 

t 

* 

I 

1 521 .7381 (1 40.4432s) 

s 2 40.6222s 40.2601 

100.6261 (1 40.4118 s) 
s 2 +0. 6248 s +0.1 936 

II 

1506*11C4(1 40.6605s) 

99.8512 

p 

s '40.6624s 40.23 C4 

s 40 0 081 8 

III 

2009.8716 

s 40.0809 

98.731 2(1 +0.31 21 s) 
s 2 40. 621 2s 40. 2025 

IV 

L 

2007. 081 3 

S'KD.0826 

99.2213 

: s4-0 . 0832 
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Table 4.1 (c) 



Turbi ne Torq ue 




















Table 4- .2 

Steam Flow Rate Through Nozzle Chest 
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XI 


5-3287 

s + 0.08+0 
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Table if .3 

Settling time (II order trans. func.) end Time 
Constants (I order trans. Tine*) 


\ Hespon -0 ] 

se 5 R.T. 1 

Model — -4 i 

1 

L.P.T.T. J 

H.P.T.T.J 

I 

S.F.R. 


11 






order 

1 1 J+ 8 

12.85 

11.98 

11.50 







JL 

II 


1,10 — 




order 

n *71 

1 2.80 

12.65 



II 


- — 



II 

order 

12.08 

12.1 7 

1 2.10 

11.90 


I 




: - 


order 

1 1 .99 

1 2.22 

12.20 



mm 


1 1 


* 

in 


12.23 

1 2.36 

12.18 



11 




- 


order 

1 2.40 

12.88 

B 



I 




* 

IV 

order 

12.20 

12.10 

12.1 9 



I 




- 


order 

1 2.51 

1 2.02 

12.30 



Abbreviations : R.T. - Reheater Temperature 

L.P.T.T.- Low Pressure Turbine Torque 
H.P.T.T.- High Pressure Turbine Torque 
S.F.R. - Steam Flow Rate# 
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Con put er Proemsing Unit lime (sees.) 


Input 

^\variablo 

T f ' ~ "5 'I 

R.T. 5 L.P.T.T. I H.P.T.T. 5 S.F.R. 

15 5 

I 5 5 

I 

92.1 3 

91 .18 

92.22 

60*38 

IX 

55*25 

55.51 

55. 61 

36 .if 3 

III 

55-90 

1 

55.72 

1 

55.11 

- 

IV 

! 

35-28 

! 

35.09 

35.11 



Abbr oration s : R.T. 


L.P.T.I. 

II.P.T.T. 


Reheat or Temperature 
Low Pressure Turbine Torque 
Pligh Pressure Turbine Torque 
Steam Plow Rate. 


S.F.R. 
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OL4PTIH V 


C ONCLUSION f LIMITATIONS OF THE PRESENT METHOD MD 
SCOP 3 FOR FUTURE WQBK 


5«1 Conclusion 

In the present work, a method has been developed 
to obtain low order models of large multi' variable systems. 
The method is applicable to both linear as well as non- 
linear systems aid does not require the transfer function 
or vector differential equation of the system to be 
specified. Only measured input- output data of the system 
at discrete intervals of time is required to determine 
its linear, low order, discrete- time model. The input- 
output data is put in the form of a svstem of difference 
equations, with model parameters as unknowns. The equa- 
tions are solved for these unknowns by the use of matrix 
pseud o- in ver se • The parameters thus determined, minimize 
the sum of the squares of the errors between the responses 
of the original system and Its low order model at the 
sampling instants. A recursive algorithm is used for 
the solution of the system of difference equations by 
matrix pseudo-inverse technique. This reduces the 
computational and data storage requirements and enables 
successively refined model of the system to be determined 
as new data from the system becomes available. 



72 


An iterative technique developed by Smith £.591 
is used to obtain the corresponding continuous- time 
model of the discrete- time model obtained by the above 
method. 

The method was applied to a nuclear reactor 
turbine system £623, £633* which has 11 non-linear 
differential equations. Reheater temperature, low 
pressure turbine torque and high pressure turbine torque 
responses for 3 * 3 % step change in main throttle valve 
area and 10 $ step charge in by-pass valve area and steam 
flow rate through nozzle for 3 . 3 % step charge in main 
throttle valve area were used to obtain the low order 
models for the system for these responses. The different 
models, with both transfer functions of II order, both 
transfer function of I order as well as combination of II 
and I order transfer functions were obtained. Responses 
of these models have been discussed and particular models 
most suitable for on-line and off-line applications have 
been suggested. 

5*2 Limitations of the present method 

The limitations of this method are as follows, 

( 1 ) Steady state error between the model and system 
responses cannot be predicted before obtaining the model. 
This is because the method minimizes the sum of squares 
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of the errors between the system and model responses 
and " fixed expression for steady state error, like 
the one given in Davison’s method^ 25 J cannot be 
derived > 

( 2 ) Least squares criterion may not giveifc the best 
results in all systems* 

(3) Values obtained from the iterative algorithm may 
not converge easily. 

0+) The time constants of all the transfer functions 
may not be the same. Thus order of the low order model 
may not be what was specified. 

5»3 Scope for future wor k 

Several suggestions are given below which may 
make the present method more useful s 

(1) In the present work, the method was applied to 

a nuclear reactor turbine system. It should be used to 
obtain low order models of several other systems to give 
more confidence on its applicability. 

(2) The method may be extended to be able to deal with 
systems -whose parameters vary slowly with time. 



(3) A t package’ computer program can "bo developed xdiich 
should directly give low order , continuous time model 
"by doing the low order modelling and converting the 
discrete- time model into the corresponding continuous- 
time model. This program should also select the most 
appropriate low order model for a given application 
by considering the response errors and model generation 
time . 

0 +) The method may be modified for 'control law’ 
reduction applications. 
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APPENDIX 


HAffilX P IElDQ-IHVERSE 


1 » Concept of h'atrix Pseuo Inverse 

For a rectangular matrix, an inverse does not 
exisb in the ordinaiy sense. It is impossible to have 
a unique and exact solution for a set of simultaneous 
equations when the number of equations are not equal to 
the number of untoowns. Two cases might arise : 


Case (a) s Number of equations exceeds the number of 
un to own s • C cn sid er the e quati on , 


H x = y (n > m) 

w* 

(nxm) (mxl ) (nxl ) 

where, 

H is a (nxm) matrix 
s is a (mxl ) vector 
and j is a (nxl ) vector , 


(JU1) 


H = 

- r 

t 7 = 

' 10' 


- 2 . 

•a* 

4- 


( A-2) 


then eqn.Ou-2) implies either 
x = 1 0, or x = 2 , 
and we cannot find a solution for x. 



The problem, then, is to find a (mxl ) vector 
x, which approximately satisfies eqn.(A-1)* The solu- 
tion is given by, 


S’ = H 31 Z (A-3) 

L 

H is a non-vnique (mxn) matrix called Left 
Inverse of H such that, 

H L H = I m (A-if) 

where, 

I m is an (mxm) identity matrix* 
letting H-k = j[ h^ h^,"! for the system of eqns* (2) we 
find that, 

h -j +ai 2 = 1 (A-5) 

Therefore H L is given by all points cn the line H-H in 
Fig. A-1 ; hence H L is not unique* The vector His 
shown lying along the line P. In other words, an infinite 
number of x' s exist, none of which are solutions of 
the original problem, 
x* = H L £ 

However, there exists a unique solution for the minimum 
LM 

Left inter se, H , given as, 

A LM 

2 = H y (JL.6) 


jfM 


exists when H has rank m or 


H 


■LM _ 


= (H H) 


-1 T 
f H 


greater, and is given by, 

(A-7) 
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We will show that the distance between Ex and 
y is minimized for the !t bost fl solution x s 


II BX - yll > ij H x - 7j| 

As an example , take systems of eqn. (A-2), 


(A.8) 


here , 


T 2 2 

H H = (1 ) + (2) = 5 


h LM 

A 

X 


= (0.2) Cl 

= Co. 2 




H x for this case is given by point in Fig. A-1 . 

P Q among all points on pp, is closest to the tip of 
h lm lies at the intersection of the lines HH and PP. 

To prove Inequality (A-8) , we let x be any 

y\ 

solution and x be the minimal norm solution, and start 
with the equality, 


i|Hx - Zji = j| H (x - x) + Hx - y|| 

= l|H (x - x) + (HH^ 1 - 2 n ) Zll 
= ||H (x - x) + (HH LM ~ f; n ) y 1| + 

2[H(£-S)i r (EH LM - ; : n )y (A.?) 

The last term of right-hand side of eqn, (A-9) vanishes; 
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T T 

£h(x - x)] (HE LM - !„) 2 = (X - X) H (HE 1 " - T n )Z 


A T I LM T v 
= (x - x) (H HH - H )2 

A. T r- T 1 TV 

= (s - x)- Lh- - 


= Oo* 

Therefore, eqn. (A- 9) reduces to, 

* T M 

ilffic - zj| =l|H(x - x) M + j{HH y-^jj 

= It H (x - x)|| + || H x - Z\\ (A-10) 

Since the first term cn t he right-hand side of this 
equation is non-negative, the inequality (A-8) holds. 

Case (b) The number of unknowns, exceeds the number 
of in owns* Consider the system of equations , 

H x = 21 (n>m) (A-11) 

(mxn) (nxl ) (mxl ) 

We cannot solve for x, because, the determinant 
t H | = 0* The problem, then, is to find a (nxl) vector 

x such that eqn* (A-11) is satisfied* 

The solution to this problem is given by, 
x = H* z ( A-12) 

The (nxm) matrix I? in this expression is called the 
■night Inverse of H* H exists when the rank of H is m 
or greater and is determined by, 
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hh r - 


1 

-in 


where. I is the (mxrrh identity matrix* 

'“ill 

TD 

As in the case (a), H" and hence the solution 
(A-12) is not unique. There exists, however, a unique 
solution, for which the norm of x is minimal. Such 
a solution x , called the Minimal Solution, is determined 

■n - r 

by the Minimum Right Inverse, H^ 1 ' 1 so that, 

A. OM 

X=K z (A.13) 


where , 

H® = H 1 < H H 1 )" 1 






